Fix trace-water phase split causing pump PS-flash failures - #4320

Review of c4e6d0b335d4a6e01a1194379d0f90146c477834, including the relationship to #4306 and the flash/performance roadmaps.

Recommendation: retain this focused repair, but request stronger public-path and performance qualification before merging. It repairs a real non-conservative endpoint without weakening PS tolerances. It is not yet evidence that the upstream phase-lifecycle cause is fixed, or that all flash workloads retain their speed.

What I verified

The six workflow groups for this head are successful. I inspected the build jobs: Java 8/21 on Ubuntu/Windows, four slow-test shards, model specifications and Javadocs passed (run).

I also compiled the exact base/head TPflash source with Java-8 language compatibility and ran isolated diagnostics on Java 17.0.20 against the same compiled NeqSim core. The base source blob matches GitHub; all 3,534 production Java files in the supporting source JAR match the local source, and connector comparison confirms that TPflash is the only production Java difference between that source and this head. These are focused diagnostics, not an independent full Maven build or reproduction of the private plant snapshot.

For the public synthetic heptane/water case at 307 K, 10.5 bara and 0.003 water mole fraction, after imposing the test's corrupted betas:

On the repaired states, maximum log-fugacity residual is below 1.31e-11; entropy and Gibbs energy return to the independently flashed reference within about 3.7e-12 J/K and 7.9e-11 J respectively. The 50/50 split is correctly used as a numerical seed. solveBeta() also updates compositions and fugacities; this is not simply overwriting beta and accepting it.

Corrections to the automated review

1. The reported NaN bypass is not demonstrated by this implementation. maximumComponentMaterialBalanceResidual() (lines 3882–3903) explicitly returns positive infinity for non-finite feed fractions, betas, compositions or accumulated inventory. Infinity passes the repair trigger and fails the acceptance comparison. Injecting NaN into the feed produced residual=Infinity, accepted=false. An explicit finite check could clarify the contract, but should not be described as fixing an established NaN bypass here. Keep the distinction from [AI-Assisted] fix(flash): preserve non-finite stability recovery #4174: a recoverable non-finite intermediate stability-monitor signal is different from accepting a non-finite final equilibrium.

2. The reflective test does exercise the new trigger. It calls rescueWaterRichEndpoint(), which calls shouldRefineTraceWaterAqueousEndpoint(). I removed only the new trigger in a diagnostic build, leaving the repair loop intact: the SRK residual stayed at 0.996999600227. Thus the statement that deleting the trigger would leave this regression passing is incorrect. The valid remaining concern is coverage of the complete public cleanup/return sequence.

Changes/qualification requested for this PR

1. Add a genuinely failing public TP → PS/pump reproducer. Keep the useful private regression. Merely replacing its reflective invocation with ThermodynamicOperations.TPflash() is insufficient: my public rerun of the same corrupted-beta case already passes on the base for SRK, PR and CPA, with component residuals around 2e-14 and unchanged entropy/Gibbs energy. Reduce the original failure to a public/synthetic fluid and state sequence that fails without this fix. Include the small-temperature scan, cold/reused operation paths, pressure change, entropy closure and pump mass/energy closure. Prove the test fails when this repair is removed. The private snapshot replay reported in the description is useful author evidence, but cannot substitute for a reproducible repository regression.

2. Distinguish a conservative phase split from a stable phase set. The new branch (lines 2285–2298) repairs the existing two phases and returns before stabilityAnalysis(). When the reference is non-conservative, shouldAcceptWaterRichCandidate() intentionally skips comparing its Gibbs energy and accepts a balanced equilibrium. Avoid comparing Gibbs energies at different inventories—that bypass is appropriate—but document that balance plus fugacity equality does not prove no third phase/lower root exists. Qualify nearby bubble/dew/three-phase points against an independent cold/stability path. If necessary, revalidate stability only after this rare invalid-endpoint repair, not on every ordinary flash. This is a qualification risk, not a claim that the supplied binary case is metastable.

3. Bound actual recovery work and make exhausted repair observable. Twenty calls to solveBeta() are not twenty Newton steps: each call can perform up to 49 internal updates, with EOS initialization inside the loop. Its normal return already follows system.init(1), so the outer candidate.init(1) appears redundant on successful iterations; verify failure paths before removing it. Record actual beta iterations, EOS init calls, clones and recovery counts. Stop on validated convergence, detect stalled/failed progress, and provide a bounded fallback or explicit failure when a known non-conservative result cannot be repaired. The current recovery caller can otherwise retain the invalid incumbent after rejection. This is an inherited acceptance/diagnostic gap exposed by the expanded repair scope.

4. Clarify tolerances and documentation. The reused production material-balance threshold is 1e-8, whereas the new test asserts 1e-10 and #2937 targets ordinary component balance at least that tight. Passing one test does not make the production gate a 1e-10 guarantee. Qualify and document the intended case-specific absolute/relative closure; do not globally tighten every solver iteration. Update the private Javadocs: they currently describe a gas/oil candidate, phase appearance/removal and a return requiring removal of one phase, which no longer describes this new branch. Add the documentation-impact assessment and a short update to the flash algorithm guide.

Speed evidence and required gate

The normal new cost for an already balanced trace-water aqueous endpoint is a component-balance scan; it does not launch the cloned repair. Preserve that property.

I ran three alternating base/head fresh-JVM pairs for complete repeated heptane/water flashes: Java 17.0.20, 256 MiB fixed heap, SerialGC, -Xbatch, 1,000 warmups, five batches of 256 flashes per case. Per-fork median thread CPU, microseconds/flash:

SRK allocation was unchanged at about 65,000 bytes/flash; PR was about 64,648 → 64,568. Post-batch balance/fugacity diagnostics stayed equal at the recorded precision. These small shared-runner diagnostics establish neither universal slowdown nor universal speedup. CPA recovery was checked numerically; CPA timing was not qualified.

Before declaring speed maintained, use paired full-flash and representative pump/process benchmarks covering dry gas, gas/oil, undersaturated trace water, established oil/aqueous, three phases, CPA and changed-state/recycle workloads. Report cold/warm/nearby-state time, allocation, recovery frequency and tail latency. Separate unchanged valid workloads from previously invalid outputs: an incorrect “fast” flash is not an acceptable performance baseline. Keep timing out of ordinary pass/fail unit tests.

Relationship to other flash work

Important current dependency status: #4306 head ec8611f is open and its build workflow is failing, despite its earlier local validation. The fast job reports sour-gas consistency, tiny-aqueous seed rejection and sulfur-filter component-balance failures. The slow job reports saturation pressure 105.892988 → 41.005951 bara. Investigate these with balance, fugacity, Gibbs/TPD and model evidence; do not simply relax assertions or change phase expectations. Do not make #4320 depend on an unqualified merge of #4306.

Recommended implementation sequence

1. Finish this bounded repair with the public failing regression, accurate documentation and no-regression performance evidence. It can remain a separate fix.
2. Resolve and qualify [AI-Assisted] refactor(flash): general constrained equilibrium and owned state #4306's lifecycle and screening changes, then test the combined source. Measure whether terminal recovery frequency falls; retain a bounded safeguard for residual hard cases.
3. Under [Campaign] TP-flash stability, multiphase equilibrium and performance roadmap #2937, progressively consolidate the repeated terminal rescue/refinement chain into one explicit contract: conserved feed and valid phase storage → warm-start solve → residual checks → stability/phase-set update when needed → bounded recovery → accepted state or explicit failure. Do this incrementally with frozen regressions, not as an unmeasured rewrite.
4. Improve speed by eliminating repeated failed trials, duplicate initialization and unnecessary cloning. Use feasible safeguarded phase-fraction solves and switch from robust successive substitution to Newton when appropriate; reuse warm-state/stability information only with clear invalidation on T/P/composition, EOS, mixing-rule and topology changes.

The numerical direction is consistent with Michelsen's stability/phase-split framework (Part I, Part II) and feasible-region safeguarded multiphase Rachford–Rice work (Okuno et al., 2010). Their inner-solver guarantees must not be mistaken for a global stability guarantee for the complete EOS flash.

The goal is fewer invalid states and fewer recovery attempts, while keeping cheap, well-qualified paths cheap—not ever more unconditional stability calculations or looser thermodynamic tolerances.

Implemented the review follow-up in 645e5e6. The original recovery remains narrowly scoped; ordinary TP/PS tolerances are unchanged.

The repair now:

• excludes forced-phase, phase-locked and solid-check states from the new existing-aqueous path;
• uses an explicit budget of 20 complete beta solves, rejects non-finite solve corrections, and removes the redundant initialization after solveBeta();
• requires local material balance 1e-10, beta/composition normalization 1e-12, and log-fugacity closure 1e-8, alongside the existing candidate-selection checks;
• checks the final endpoint after an attempted repair, allowing subsequent phase topology changes, and throws an informative exception if it remains invalid.

The final check is local numerical qualification, not a global stability certificate. Balanced endpoints incur no extra candidate clone or EOS initialization in the guard.

Validation

130/130 focused tests pass, including public SRK/PR/CPA temperature sweeps, pressure changes and return, reused operations, cold-reference agreement, and actual pump entropy/mass/energy closure. Injected-state tests cover strict acceptance, failure reporting and the cheap guard. The original corrupted-endpoint regression remains.

The same 130 tests pass in an isolated combined build with the targeted #4306 lifecycle fixes at 518867a. Those fixes preserve root selectors and physical phase mappings during rollback, index Newton beta logically, and reject invalid higher-Gibbs aqueous seeds. #4306 has its own 145-test passing focused suite, but its 3 sour, 3 wet-sulfur and 1 scrubber failures remain explicit merge blockers. Slower sour-polishing and broader Newton-objective experiments were excluded.

Java 8-language compilation used ECJ; execution used Java 17 and verified baseline compiled dependencies. This was not a full local Maven/JDK build. Direct Spotless apply/check, documentation search, engineering inventory generation, diff checks, and both pre-commit/pre-push stages pass. The new GitHub build/test/Javadoc matrix remains authoritative.

Performance evidence

Baseline: this PR's previous exact head c4e6d0b; candidate: 645e5e6. Identical harness, Linux Java 17, -Xms256m -Xmx256m -XX:+UseSerialGC -Xbatch, three alternating fresh-JVM pairs.

The complete-flash screen covers 12 families × unchanged/changed states, 5 batches × 32 measured flashes per state, 64 warmups (4096 for the three cheap families). All 720 batch-end numerical snapshots pass; all 360 matched baseline/candidate snapshots are identical. Median allocation differences are at most 4 bytes/flash. An injected SRK repair also produced exactly the same balance/fugacity residuals while reducing counted init(int) calls from 43 to 42.

Short-warmup timing was noisy and showed increases in some CPA/condensate cases. Longer-warmup checks were therefore run on a fixed CPU with 1000 warmups, 5 measured batches, and three alternating JVM pairs:

The trace-water repeated-equilibrium checks used 1000 warmups and the same JVM flags, but were not CPU-pinned. Individual paired measurements still varied considerably; these are representative measurements, not a universal throughput or tail-latency guarantee. The longer-warmup checks did not reproduce the screen's CPA/condensate slowdown.

The existing large-process harness also passed all 10 measured serial and 10 recycle changed-state samples on both versions (size 4, sequential strategy, 10 warmups). Single-pair median elapsed times were serial 117.0 → 92.2 ms and recycle 311.2 → 318.8 ms; these process timings are diagnostic only.

The algorithm guide and task log document the trigger, exclusions, acceptance contract, retry cost and validation limits. The original private process snapshot was unavailable during this follow-up; public cases qualify lifecycle behavior and do not establish full-plant convergence. Neither PR has been merged.

Updated this PR in be5d85f, then merged current master (b0de5165) in 2b67ee0 to clear the conflict. GitHub now reports no merge conflicts.

The CI failure was caused by subtractive cancellation in SystemThermo.setTotalFlowRate: reservoir-sized component inventories were subtracted down to a small stream, corrupting the overall composition. The public reservoir regression reproduces the exact old-head repair exception. Rebuilding positive normalized component inventories fixes the feed itself, preserving trace species.

The performance changes retain the qualified gas/aqueous candidate instead of cloning/reflashing it, and stop non-finite or repeatedly stalled repair iterations early. The strict normalization (1e-12), material (1e-10) and log-fugacity (1e-8) acceptance gates, phase exclusions and final failure reporting remain intact. Ordinary TP/PS tolerances are unchanged.

Validation

• Exact old head 645e5e6: the two regression classes expose seven failures, including the public reservoir path and the duplicate-clone assertion.
• Final merged source: 530 tests, 0 failures, 0 errors, 2 existing disabled tests (528 passed). All previously failing reservoir, well and field-development cases pass, including the slow lifecycle cases. Added master-side network and wax regressions also pass.
• Full Maven production/test compilation on OpenJDK 17.0.20; changed Java files also compile with --release 8 (language/API compilation, not a Java 8 runtime claim).
• Direct Spotless apply/check, documentation audit, engineering-coverage check, pre-commit, pre-push, whitespace checks and Javadoc generation all exit 0. Javadoc needed a runtime-only launcher for the installed JDK module because the executable is absent in this environment.
• Both existing review threads are addressed and resolved.

Commands, run after the NeqSim Maven bootstrap:

./mvnw -B -ntp test '-Dtest=TPflash*Test,TPmultiflash*Test,PSFlash*Test,PHflash*Test,SystemThermo*Test,PumpTest,PumpSinglePhaseLiquidPSflashTest,ReservoirSurveillanceTest,WellFlowTest,FieldDevelopmentNPVTest,NorwegianOilFieldLifecycleTest,NamedComponentWaxFractionTest,WaxFractionSimTest,FieldNetworkTopologyTest,FieldWellNetworkProcessUnitTest,LoopedPipeNetworkTest' -DexcludedTestGroups= -Djacoco.skip=true
python3 devtools/run_spotless.py apply
python3 devtools/run_spotless.py check
python3 devtools/check_documentation_search.py
python3 devtools/build_engineering_coverage.py --check
pre-commit run --all-files --hook-stage pre-commit
pre-commit run --all-files --hook-stage pre-push
./mvnw -B -ntp javadoc:javadoc -DskipTests

Performance and numerical evidence
To isolate this fix, paired benchmarks compare exact old production classes 645e5e6 with be5d85f; final master integration was validated separately by the 530-test run. Three fresh-JVM pairs, alternating version order, each pair pinned to a CPU; later pairs ran concurrently on separate CPUs. OpenJDK 17.0.20, -Xms256m -Xmx256m -XX:+UseSerialGC -Xbatch.

The injected diagnostic explicitly sets the equilibrated methane/water gas/aqueous betas to 0.5/0.5 before timing the private repair (1,024 warmups, seven batches of 256 repairs; 42/42 measured endpoint snapshots pass). It measures the removed clone/reflash; it is not a complete public-flash speedup claim. The counted regression independently proves one candidate clone and one live initialization.

The public TP harness covers 12 families × unchanged/changed conditions × 5 batches × 3 pairs, with 32 flashes/batch and 128 warmups (4,096 for the cheap cases). Maximum material residual is 4.68e-12, log-fugacity residual 5.24e-11, phase normalization 3.34e-16, beta normalization zero. Per-family median CPU changes: dry gas +1.1%, hot gas −4.0%, two-phase −1.0%, CPA water +2.6%, condensate +0.6%, cold condensate +1.0%, multiphase condensate −0.3%, trace aqueous SRK +0.3%, PR −1.4%, CPA +1.0%, undersaturated +2.4%, three-phase +4.6%. Individual short-batch outliers occur, so no general flash speedup is claimed.

The initial short-warmup recycle result was noisy (+6.4% median, −1.0% to +26.0%); the longer-warmup check above resolves that uncertainty. Process iteration counts are unchanged, and the harness's existing mass/component/energy criteria all pass. This does not tighten the process recycle tolerance or claim private-plant convergence. Invalid old reservoir results are never counted as a fast baseline.

Reproduce the public harnesses with each version's compiled classpath:

java -Xms256m -Xmx256m -XX:+UseSerialGC -Xbatch -cp "$CP" neqsim.thermodynamicoperations.flashops.TPflashEndpointRecoveryBenchmark /tmp/flash.json 32 5 128
java -Xms256m -Xmx256m -XX:+UseSerialGC -Xbatch -cp "$CP" neqsim.process.processmodel.LargeProcessSteadyStateBenchmark serial changed 10 10 4 sequential /tmp/serial.json
java -Xms256m -Xmx256m -XX:+UseSerialGC -Xbatch -cp "$CP" neqsim.process.processmodel.LargeProcessSteadyStateBenchmark recycle changed 100 30 4 sequential /tmp/recycle.json

CI: Exact-head run 38539212902. Model-specification jobs and Java 21 Javadocs have passed; full Ubuntu/Windows and slow-test jobs are still running. Formatting, pre-commit, documentation search, engineering coverage, CodeQL security and code-quality checks have passed. This remains an open PR; it has not been merged.

Published the CI regression repair at 15d1b8defe7cff39899ff8592f8ed75baebf5e2f, using an expected-head lease against 2b67ee03629bf845687d14e8a985be982fd4f8a0.

The previous head reproduced eight assertion failures and two errors. An independent exact-master run passed all 76 tests in those seven classes, confirming these were branch regressions. The common cause was applying positive inventory reconstruction to every flow update; ordinary floating-point endpoints and empty-stream behavior changed. The repair limits reconstruction to downward updates whose component or total difference update would lose more than 1e-12 relative accuracy. Upward, same-rate and accurate ordinary updates retain their established arithmetic. Extreme 1e-12–1e18 mol scaling and the 1e-10 trace fraction remain tested.

The two-fluid VLP failure had a separate direct cause: scaling its property-only fluid clone to a near-zero trial rate. That scaling is unnecessary because the calculation reads intensive properties and already applies the trial rate to its velocities. Removing the setter call restores positive near-zero continuity and avoids one init(0) plus one init(1) per positive evaluation. New assertions check reference inventories from 100 to 1e6 Sm3/day and verify the source stream is unchanged. The explicit zero-flow boundary is retained.

The reservoir reference change is supported by a material audit, not a widened tolerance. The old oil producer was 0.0322% above its requested rate (1,260,096.07887 versus 1,259,689.92248 kg/hr). Correcting that rate changes the ten-year pressure to 355.175066 bara. The existing 0.1 bar tolerance is retained, and the test now asserts all four stream rates, component inventory against integrated production/injection, and constant volume at every annual step. Independent maximum relative component imbalance was 1.50e-16; relative volume error was below 1.5e-12.

Validation on the published source tree

• Expanded Maven run: 607 reported = 605 passed + two pre-existing disabled; zero failures/errors, slow exclusions disabled, build successful in 6m35s.
• Original failing classes all pass: adsorption 23, distillation 26, exchanger 4, reservoir 1, well 15, constraint activity 3, operating action 4.
• Existing/new inventory cases: 6 pass. Broader coverage includes TP/PS/PH, cubic/CPA/lifecycle/aqueous endpoint cases, pumps, networks, reservoirs and flow units.
• All five edited Java files compile with --release 8. Maven/test runtime is OpenJDK 17.0.20, not a Java 8 runtime certification.
• Direct Spotless apply/check; documentation search (788 Markdown + one HTML); engineering inventory generation/check; pre-commit and pre-push hooks; whitespace checks; Javadocs: all pass. Javadocs use the installed JDK module through a runtime-only launcher because the environment has no standalone javadoc executable. No hook bypass.

A focused reproduction command for the affected and inventory suites is:

./mvnw -Dtest=AdsorptionBedTest,DistillationColumnTest,MultiStreamHeatExchangerTest,SimpleReservoirTest,WellSystemTest,ConstraintActivityAnalyzerTest,ProcessModelOperatingActionTest,SystemThermoFlowRateInventoryTest -DexcludedTestGroups= -Djacoco.skip=true test

Performance limits: three fresh-JVM setter pairs and three warmed 10,000-iteration process pairs had substantial timing variability. The warmed process CPU differences ranged from −20.3% to +31.1%; allocation varied by less than 0.05%, and all 60,000 measured process iterations were valid. These measurements do not establish a speedup or a regression. The defensible optimization is removal of the two unnecessary VLP initialization calls; the cancellation guard adds no EOS initialization, clone or temporary array. The earlier trace-water candidate-transfer optimization is unchanged.

Related clone/root-selector problems are being repaired independently in #4306; its changes are not copied into this branch. GitHub CI on this new exact head is pending. No merge performed.